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ABSTRACT 

NGC 4494 is one of several intermediate-luminosity elliptical galaxies inferred to 
have an unusually diffuse dark matter halo. We use the x 2 -made-to-measure particle 
code NMAGIC to construct axisymmetric models of NGC 4494 from photometric and 
various kinematic data. The extended kinematics include light spectra in multiple 
slitlets out to 3.5i? c , and hundreds of planetary nebulae velocities out to ~ 7R C , thus 
allowing us to probe the dark matter content and orbital structure in the halo. 

We use Monte Carlo simulations to estimate confidence boundaries for the halo 
parameters, given our data and modelling set-up. We find that the true potential of 
the dark matter halo is recovered within AG(merit function)< 26 (A% 2 < 59) at 70% 
confidence level (C.L.), and within AG < 32 (Ax 2 < 70) at 90% C.L.. These numbers 
are much larger than the usually assumed A% 2 = 2.3(4.6) for 70% (90%) C.L. for two 
free parameters, perhaps case-dependent, but calling into question the general validity 
of the standard assumptions used for halo and black hole mass determinations. 

The best-fitting models for NGC 4494 have a dark matter fraction of about 0.6±0.I 
at 5i? e (70% C.L.), and are embedded in a dark matter halo with circular velocity 
~ 200kms~ . The total circular velocity curve (CVC) is approximately flat at v c = 
220kms~ outside ~ 0.5i? o . The orbital anisotropy of the stars is moderately radial. 
These results are independent of the assumed inclination of the galaxy, and edge-on 
models are preferred. Comparing with the halos of NGC 3379 and NGC 4697, whose 
velocity dispersion profiles also decrease rapidly from the center outwards, the outer 
CVCs and dark matter halos are quite similar. NGC 4494 shows a particularly high 
dark matter fraction inside ~ 3R e , and a strong concentration of baryons in the center. 

Key words: galaxies: kinematics and dynamics - galaxies: elliptical and lenticular, 
cD - galaxies: individual: NGC 4494 - galaxies: halos - cosmology: dark matter - 
methods: N-body simulations 



1 INTRODUCTION 

The formation and evolution of elliptical galaxies has been 
an important issue in extragalactic astrophysics since a long 
time. Being collisionless to a very good approximation, ellip- 
ticals retain relics of their formation history in the present- 
day orbital structure, especially in their halos due to the 
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longer dynamical time scales. According to the currently 
favoured hierarchical formation scenario in a ACDM cos- 
mology, these halos are dark matter dominated. The ambi- 
tious task of inferring both the orbital structure and mass 
distribution of ellipticals is commonly tackled by dynamical 
modelling of the observational data. 

Unfortunately, the lack of an ubi quitous tracer such 
as HI gas in spiral galaxies (but see [ Bertola et all 1 19931 : 
iFranx et al.l 1 19941 : lOosterloo et al.l 120021 ') makes mass mea- 
surements in elliptical galaxies challenging. The strongest 
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evidence for dark matter halos is found for giant ellipti- 
cals whose mass distribution c an be determined from X-ray 
emission of the hot g a s {e.g. Loewenstein fe White! 1 19991 ; 



Humphrey et al.ll2006l ; iDas 'et al 



2010a [) or strong gravita- 
iKeetonl 



Maoz fe R: 



Auger et al 



1993 



20101 ). These 



tiona l lensing techniques (e.g. 
l200ll ; iTreu fe Koopmans! |2004| ; 
studies are consistent with m assive dark halos, and nearly 
isothermal total mass profiles (|Ciotti et al.ll2009l ). 

By contrast, the situation with less massive, X-ray faint 
ellipticals is more controversial. However, dynamical mod- 
els, particularly when fitting higher order moments of the 
line-of-sight velocity distribution (LOSVD) from stellar ab- 
sorption lines, eventually ascertained the presence of dark 
matter halos around these intermediate-luminosity ellipti- 
cals and are generally consistent with flat circular velocity 
curves. With LOSVD measurements limited by the rapid 
fall-off of the stellar surface brightness, the kinematics of 
discrete tracers such as planetary nebulae (PNe) and glob- 
ular clusters (GCs) usually represent the only possibility 
to probe the mass distribution and orbital structure be- 



1995 




Mcndez et al.l l200ll; iDouglas et al.l |2002|; iPeng et al. 


2004 


de Lorenzi et all 2008. 20091; Na 


oolitano et al. 


2009; 


Das et al.ll2011; Naoolitano et al. 2011; 


Deason et al.l 


2012). 



Curiously, the PNe velocity dispersion profiles of some 
of the nearby intermediate-luminosity ellipticals show a 
strong, quasi-Kepler ian decline with radius outside lR e 
dCoccato et al.l I2009T), suggesting very little (if any) dark 
matter (|Romanowskv et al.l 120031 ). The aim of the present 
paper is to expand the sample of modelled quasi-Keplerian 
intermediate luminosity ellipticals, focussing on NGC 4494, 
and then compare the results with those previously ob- 
tained for other galaxi es of this c lass, n amely NGC 4697 
j2008l . hereafte r iDLOSl) and NGC 3379 
2009, hereafter [PL09l ). 



de Lorenzi et al 



de Lorenzi et al 



NGC 4494 is an E1-E2 elliptical galaxy in the outer re- 
gions of the Virgo cluster, with a smooth ligh t profile and 
an in termediate stell ar mass of about 10 11 Mq |Foster et alj 
l201ll . hereafter iFllI ). It has been variou sly described as a 
loose group member ([Forbes et aT1ll996l ) or isolated galaxy 
|Lackner fe Ostrikeill2010l) NGC 4494 is classified as a fast 
rotator by Emsellem et al.l (|201ll . Ar ~ 0.2). Among the 
peculiarities of thi s galaxy, a sharp central ring of dust 
l|Forbes et al.l Il995l; Lauer et al.l 2005 | ). a kinematically de - 



coupled core ( Bender et alj 19941 ; iKrainovic et alj 20 111) 



moderate rotation ~ 60kms out to ~ 3R e (iProctor et al 



2009), and an outward-decreasing AR-profile (|Coccato et al 



2009) have been reported. The velocity dispersion profile 



of NGC 4494 decreases rapidly from ab out 160kms in 
the c enter to about 70kms _1 at ~ 7R e (|Napolitano et all 
120091 . hereafter Ho|), hinting to a possible deficiency in dark 
matter. Also, the X-ray flux of NGC 4494 is two hundred 
times fainter than that of other galaxies of the same opti- 
cal luminosity, which has been interpreted as the result of 
a recent interaction which dep leted the gas, or again as ev- 
idence for little dark matter jO'Sullivan fe Ponmanl |2004| ; 
iFukazawa et al.ll2006l ). 

NGC 4494 has b een recently modelled using 
spherical Jeans models iloi), and axisymmetric par- 
ticle models constructed with the iterative method 
dRodionov fe Athanassoulal 120111) . The available observa- 
tional data, which included PNe velocities out to ~ 7R e 



(Hoi), were best fit by low concentration dark halos, 
with some uncertainties related to the adopted modelling 
assumpt ions, the differences between the fits for different 
models (IN09T), and the limited n umber of explored models 
|Rodionov fc Athanassoulall201ll) . 

Recently, new observational data consisting of stellar 
absorption line kinematics in mul tiple slitlets out to ~ 3.5R e 
became available for this galaxy (|Proctor et al. IHoqIIfIi]). 
Moreover, in the work of iDeason et al.l (|2012l ) NGC 4494 
appeared as an outlier with curiously low dark matter frac- 
tion within 5R a , with respect to model predictions assuming 
either a Salpeter or a Chabrier initial mass function. These 
facts prompted us to undertake a further careful analysis of 
the dark matter content and orbital structure of NGC 4494, 
incorporating as many observational data as currently avail- 
able and assessing the uncertainties in the recovery of dark 
halo parameters via dynamical models. 

In this paper, we construct new dynamical mod- 
els fitting all available photometric and kinematic data 
with the flexibl e part icle code NMAGIC (|de Lorenzi et all 
120071 , hereafter IDL07I ). which implements a modified ver- 
sio n of the made-to - measu re (M2M) technique proposed 
by ISver fe Tremainel (|l996l ). suitable for the modelling of 
observational data with errors (X 2 M2M). NMAGIC works 
by slowly correcting the particle weights of an evolving N- 
body system until a satisfactory compromise is achieved 
between the goodness-of-fit to the observational data, and 
some degree of regularization of the underlying particle 
model. M ore recen t impl ementations of the method can be 
found in iDehnenl ([2009) , who p roposed a different tech- 
nique for the weight adapta tion, lLong fe Mad (|2010l ). and 
iMorganti fe Gerhard! (|2012l . hereafter |MG12|) , who intro- 
duced a Moving Prior Regularization method to generate 
smooth x 2 M2M particle models fitting noisy data without 
erasing the global phase-space structures. 

So far, the M2M technique has been used to inves- 
tigate the dynamics of the Milky W ay's bulge and disk 
|Bissantz et al.ll2004l ; lLong et alj|2013h . the mass distribu- 
tion and orbital structure in the oute r halos of elliptical 
galaxies (|DL08l ; |PL09l ; iDas et al.lboill ). and the dynamics 
of a sample of SAU RON elliptical and lenticular galaxies 
l|Long fe Mad |2012| ). In this paper, we construct axisym- 
metric NMAGIC models for NGC 4494 which include PNe 
(N09) and new stellar absorption line kinematics in slitlets 
l|Flll) . for different dark matter contributions to the total 
gravitational potential, and for different inclinations. 

One of the key points in this work is to study the accu- 
racy with which the parameters of our models can be esti- 
mated given the observational data at hand. M2M particle 
methods work by adjusting the weights of a large number of 
particles in order to achieve a good match to the data. Sim- 
ilarly, in the more common Schwarzschild methods a large 
number of orbital weights is adjusted to fit the data. Typ- 
ically, there are many more weights than data constraints, 
hence the number of degrees of freedom^ is likely to be much 
smaller than the number of data points (|Cretton et al I I2OOOI ; 



1 The number of degrees of freedom is defined as the number 
of constraints (data points plus constraints introduced by e.g. 
regularization) minus the number of free parameters (model pa- 
rameters plus fitted weights). 
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iGebhardt et~ai1l2000l ). However, to determine the effective 
number of free parameters involved in the modelling is very 
difficult. 

It is common practise to use the relative differences 
Ay 2 = X 2 ~ Xmin to measure confidence limits for the 
subset of parameters specifying the mass model. On the 
assumption that th e observational errors are Gaussian 
distributed (but see Ivan der Marel etTail l200Ch . that the 
model is linear in the parameters, and the number of 
degre es of freedom is positive, A% 2 follows x 2-s t a tistics 
(e.g. iPress et al.lll992T l. Many dynamical modelling studies 
have therefore estimated, e.g., 68% confidence limits for 1 
(2,3) free model parameter(s) from A y 2 = 1(2.3,3.5) con- 
tours i n these parameters (e.g . van der Marel et al 

19991; iBarth et"al . 

20021: 



199S ; Cretton fc van den Boschl 



200 it IVerolme et all 120021; ICappellari et al 



Gebhardt et al.l 120031: IValluri et all 12004 IShapiro et all 



20061: IChaname et all 120081: ICappellari et al 
van den Bosch fc de Zeeuwl l20ld; Murphv et al. 



2009 



2011 



Adams et al.l I2012T ). but see Ivan den Bosch fc van de Venl 



(2009). Here, prompted by large Ay 2 values between quite 
similar models, we use a Monte Carlo simulation approach 
to estimate confidence intervals for the paramet ers based on 
models of NGC 4494-like mock galaxies (see als o IPress et alj 
ll992l : lGerhard et al.lll998l : iThomas et al.ll2005T ). 

The paper is organized as follows. In Section [2] we de- 
scribe the observational data used for NGC 4494. In Sec- 
tion |3] we outline the modelling technique. In Section [3] we 
construct NGC 4494-like target galaxies and their observ- 
ables in different dark matter halos, and use them to cal- 
ibrate the optimal amount of regularization, and to esti- 
mate the confidence levels for parameter estimation with 
our NMAGIC models. Dynamical models of NGC 4494 for 
a range of dark matter halo potentials and inclinations are 
then constructed in Section [5] The main implications of our 
findings are discussed in Section [6] Finally, the paper closes 
with our conclusions in Section [7] 



2 OBSERVATIONAL DATA 

In this Section we describe the observational data that we 
will use for modelling the elli ptical galaxy NGC 4494. We 
adopt a distance of 15.8 Mpc jTonrv et al.l bOQl). so that 1 
kpc=13", a systemic velo city v = 1344kms~ 1 , from NED, 
and Re = 49" » 3.77 kpc (|de Vaucouleurs et al.lfl99ll ). 



2.1 Photometric data and deprojection 

As photometric data, we use the V-band surface-brigh tness 
profile, ellipticity and shape parameter 04 values of lN09l . 
The photometric data extend to 273" along the major axis, 
and are a c ombination of HST data in the V and I bands 
inside 4.3" (|Lauer et al. 20051). ground-ba sed observations 
in BVI out to 32" l|Goudfrooii et "all 119941 ). and Megacam 
data from the Sloan Digital Sky Survey g' filter l|N09l) . The 
total exti nction -corrected luminosity in the V-band is 2.6 X 
10 10 Lv, Q (Hoi). 

A Sersic fit to the surface brightness profile outside the 
central dust region (R > 5.6") gives n = 3.30 ijXOflh. The 
observed ellipticity is e = 0.15 — 0.20 (axis ratio q = 0.85 — 
0.80) for R < Re/2, outside of which NGC 4494 becomes 




Figure 1. From top to bottom: V-band surface-brightness, ellip- 
ticity, and fourth-order isophotal shape coefficient 04 as a func- 
tion of major-axis projected radius for NGC 4494. Circles are the 
measurements, lines represent the reprojected values from ax- 
isymmetric deprojections obtained with an assumed inclination 
of 90° (black), 70°(green), and 45°(red). 



rounder, with e = 0.13 - 0.15 (q = 0.87 - 0.85) for 1 - 
1.5-Rc. The radial profiles of surface-brightness, ellipticity, 
and shape parameter 04 are shown in Fig. [1] 

In our NMAGIC models for the galaxy NGC 4494, we 
will not use the 2D surface brightness, but rather its 3D 
deprojected luminosity density. The deprojection of the sur- 
face brightness is unique only for spherical or edge -on ax- 
isymmetric systems (e.g. iGerhard fc Binnevl 1 19961 ). Here, 
we consider axisymmetric deprojections for inclinations i — 
90°, 70°, and 45°, which is very close to the minimum incli- 
nati on allowed by the observed flattening of NGC 4494 (lN09l : 
iFllI ). For each inclination angle, we use the maxi mum penal- 
ized li kelihood scheme and program described in lMagorrianl 
l|l999l ) to find a smooth axisymmetric density distribution 
consistent with the surface-brightness profile. The method 
favours a power-law density profile in the radial region not 
constrained by photometric data. 

The overall good agreement between the measured and 
reprojected surface-brightness, ellipticity, and shape param- 
eter 04 is shown in Fig. \T\ For i — 45° , the method appar- 
ently has some difficulties in reproducing the boxy isophotes 
at R ~ 200" accurately. 



2.2 Kinematic data 

We combine three kinematic data sets in order to achieve the 
widest possible spatial coverage and probe the mass distri- 
bution and orbital structure far out in the halo of NGC 4494. 
In particular, we use lo ng-slit absorption lin e kinematics ex- 
tending out to ~ 2R e (|Coccato et a l. 2009), absorption line 
kinematics in slitlets out to ~ 3. 5R e (Fill) , and PNe line-of- 
sight velocities reaching 



7R a |N09l) . The spatial coverage 
of the combined kinematic constraints can be appreciated in 
Fig.H 
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Figure 2. Spatial distribution of published kinematic data for 
NGC 4494. Long-slits along the major and minor axis are plotted, 
where the x-axis coincides with the major axis of the galaxy. The 
positions of PNe are marked by black crosses, those of the slitlets 
by green squares. Ellipses represent 3, 6, and 9R e , for an axis ratio 
q = 0.87. 



2.2.1 Stellar- absorption line slit data 

Long-slit absorption line kinem atics within ~ 2R e were pre- 
sented in ICoccato etail l|2009h . and consist of line-of-sight 
velocity v, velocity dispersion a, and higher-order Gauss- 
Hermite coefficients h-j, and Ha, along the major and minor 
axis of NGC 4494. 

The original data are shown in Fig. [3] with black dots. 
They are consistent with small or zero rotation along the 
minor axis, and substantial major axis rotation, flattening 
beyond 20" at V ~ 60kms _1 out to 120". In the inner 10" 
along the major axis, the signature of the decoupled core 
is apparent. The velocity dispersion decreases from about 
160kms _1 at the center to about 80kms _1 at ~ 100". 

We noticed a systematic offset between the velocities 
and /13 coefficients measured along the major and minor 
axes in the central arcsec. This leads us to suspect an offset 
of the minor axis slit from the galaxy centre, in the South 
direction. The required offset (less than 1") is smaller than 
the slit width and the average seeing. Furthermore, the /13 
measurements along the minor axis are overall negative on 
both sides, which is unexpected even for a triaxial system 
with minor axis rotation. For the modelling, we therefore 
replace the measurements of v and /13 along the minor axis 
with Gaussian random variates with zero mean and la dis- 
persion equal to the observational errors. 

Inside 32" there are many nearby data points whose 
respective rms deviations are larger than their error bars. 
In order to reduce their impact on the modelling, we run a 
central moving average over the data within 32", averaging 
over 7 data points (3 points on each side), and substituting 
each point with the value of the average. Because the errors 
of the neighbouring points are very similar, we decided to 
leave them at their original values. 

Then, we minimize the impact of the feature at ~ 20" 
along the major axis, which is due to contamination from a 



100 
> 
-100 

1 50 

b 

100 

50 
0.1 



-0.1 
0.2 
0.1 
^ 
-0.1 
-0.2 



40 
20 
> 
-20 
-40 

1 50 

b 

1 00 

50 
0.1 
j? 
-0.1 
0.2 
0.1 
^ 
-0.1 
-0.2 



-III 

■ ■ ■ t i ■ »«» ■ 


, , , I ' ' ' ' I ' ' ' '_ 

■■111. . , 


H 1 1 1 1 1 1 1 1 h-, 

• 1 i i i '"^ 


H 1 1 1 1 1 1 1 1 1 h 

— 

t — 1 — 1 — — 1 — 1 — 1 — 1 — — 1 — 1 — 1 — E 


H 1 1 1 1 1 1 1 1 1 — 


■f i I I I I I T 1 


H — — 1 — 1 — 1 — 1 — — 1 — 1 — 1 — 1 — — 1— 

; .... ***** 


^ 1 1 1 — 1 1 — 1 — 1 1 1 F 

$ • * J j \ 


-100 -50 


50 100 


i i i i i i i i i i i 

: : • ' • 






= i i i i i i i i i i i i i i i i i = 

I — 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 H 


- * * i » if^^kn^l i 1 } 


-100 -50 

R ( 


50 1 00 

') 



Figure 3. Long-slit kinematics along the major (top) and minor 
axis (bottom) of NGC 4494. Black and red dots represent the 
original data and the data used for the modelling, respectively. 
Sec Section |2.2. II for details. 



foreground star, by artificially increasing the error bars of 
those measurements. 

As we are interested in axisymmetric models of the 
data, we finally symmetrize the slit data set, as in iDLOgt 
In practice, we average the values of the measured kinemat- 
ics at two similar radii (R+,R-) on both sides of the slit 
with respect to the center. Since the kinematic data show 
major axis rotation, we take into account the sign reversal of 
11 and /13 when symmetrizing the major axis slit data. Then, 
the new symmetrized data point is set equal to the weighted 
mean of the points on both sides, with weights proportional 
to the inverse square of the original errors. If cr+ and <r_ are 
the errors on both sides, the new error a on the symmetrized 
data points is set equal to the maximum of 



er 2 <rl a 2 



(1) 



and half of the deviation between the original data points, 
which includes systematic errors between both sides of the 
galaxy. 

The resulting data, with their respective error bars, are 
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shown with red dots in Fig. [3] where they can be compared 
with the original measurements. 

2.2.2 Stellar-absorption line slitlets data 

We include in our kinematic constraints the measurements 
of v, a, /13, and /14 i n 11 5 galaxy light spectra in slitlets 
rece ntly published bv lFlll . which extend out to 3.5i? c - 

IF 111 discussed the generally good agreement between 
these data an d the long-slit absorption line kinematics of 
ICoccato et al.l l|200Sh . Also, they reported that their uncer- 
tainties are likely to be slightly underestimated. 

On the whole, the absorption line kinematic data do not 
show significant evidence for minor axis rotation. Therefore, 
being interested in axisymmetric models of the data, we 4- 
fold the original sample of slitlets in order to decrease the 
impact of data asymmetries in our models. To this end we 
reflect the measurement points around both the galaxy's 
major and minor axes, maintaining errors, and taking into 
account the sign reversal of v and hi,. 

Finally, we look for possible outliers in the 4-folded data 
set consisting of 460 data points. To this aim, we use the 
following procedure: for each slitlet, we compute the value 
of the average field of v, a, /13, and hi from its 20 near- 
est neighbours, excluding the slitlet itself. This average is 
worked out by excluding the lowest and highest 2 values in 
the neighbourhood, i.e. considering only the central 80% of 
the distribution. Moreover, the average is a weighted one, 
with weights equal to the inverse square of the larger of the 
individual observational error and the median of the obser- 
vational errors in the central 80% of the distribution. In the 
same way, we calculate a weighted rms a in each neighbour- 
hood. Then, we flag as outlier any point deviating by more 
than 2a from the weighted mean of its neighbours. 

With this procedure, 158 data points are flagged as 
outliers in velocity and/or velocity dispersion, and removed 
from the 4-folded data set. Of the resulting data set con- 
sisting of 302 slitlets, 68 points are marked as outliers in 
/13 or /14, and therefore we only consider their velocity and 
velocity dispersion in the modelling below. 

The slitlets data points, together with the outliers de- 
termined in this way, are shown in Fig. [4] 

2.2.3 Planetary nebulae velocities 

Our kinematic observables are co mplet ed by the 267 PNe 
line-of-sight velocities obtain ed by |N09| with the Planetary 
Nebulae Spectrograph fPNS. iDouglas et al.l 2002). As illus- 
trated in Fig. [2l PNe extend out to ~ 7 Re. 

Fig. 0] shows the comparison of all kinematic data. In 
the region of overlap, the kinem atics of PNe are consistent 
with those of the stars fsee lN09l . Fig. 6 therein). 

In order to look for possible outliers in the sam- 
ple of PNe, we use the friendless algorithm presented in 
iMerrett et all l|2003h . which flags as outlier any object de- 
viating by more than n x a from the velocity distribution 
of its M nearest neighbours, where a is the rms computed 
in each neighbourhood. Adopting n — 2.5 and M — 20, we 
remove 10 outliers from the original sample of PNe. These 
outliers are highlighted in Fig. O which shows the projected 
phase-space distribution of PNe. 




100 200 300 400 



R (") 

Figure 5. Distribution of the line-of-sight velocities of PNe as a 
function of projected radius. Red crosses mark the outliers iden- 
tified by the friendless algorithm. The PNS filter band-pass is 
overplotted in green. 

We are left with a catalogue of 257 PNe, whose size 
we double by applying point-symmetry, i.e. generating for 
every PN (x,y,v) the symmetric PN (-x,-y,-v). Such point- 
symmetric velocity fields are expected in axisymmetric (or 
triaxial) potentials. From previous experiments we know 
that increasing the number of PNe by a further factor of 
two does not improve the modelling significantly. 



3 MODELLING DATA WITH NMAGIC 

In this Section we briefly describe the NMAGIC made-to- 
measure (M2M) modelling technique. We also explain how 
the initial particle model is set up, and how the photometric 
and kinematic target observables are preprocessed for the 
dynamical modelling. 

3.1 NMAGIC 

The parallel code NMAGIC (|DL07l) is an implementation of 
a particle-based method to create x 2 M2M models in agree- 
ment with observations of galaxies. The al gorithm is a slight 
modif ication of the technique proposed bv lSver fc Tremainej 
(1996), designed to model observational data with errors. 

The basic idea behind M2M particle methods is to train 
a system of i = 1, . . . , N particles to reproduce the observ- 
ables of a target galaxy by maximizing the function 

F = -±x 2 + vS + £ (2) 

with respect to the particle weights u>i. This maximization 
strikes a compromise between the goodness of the fit (\ 2 ) 
in terms of deviations between target and particle model 
observables, and a pseudo-entropy functional 

JV 

NSH (3) 

i=i 
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Figure 4. Slitlcts kinematics in the 4-folded data set {circles), compared with long-slit {blue dots) and PNe data {orange triangles). Left 
and right column are for the measurements around the major and minor axis, respectively. Filled circles represent the data in a cone of 
±15° around the axis, while empty circles are the data points in cones 30° ± 15°, 60° ± 15°, respectively. Outliers determined by the 
method described in Section l2.2.2l are marked red. 



l|MG12t ) which serves the purpose of regularization. Finally, 
the likelihood term C is added to equation ([2]) to accoun t 
for the likelihood of a sample of PNe velocities (see IDLO&I ) . 

Since in typical applications the number of particles 
is much higher than the number of data constraints on 
the particle model, regularization is essential. In standard 
M2M, regularization is achieved by pushing the individual 
particle weights towards a smooth distribution of predeter- 
mined priors Wi, which mirror the initial particle weight 
distribution and are kep t constant during the modelling. 
As we showed in |MG12| . such a Global Weight Entropy 
scheme makes it hard to reconcile smoothness and orbital 
anisotropics in the final particle model. Therefore, in this 
work we adopt the alternati ve Movi ng Prior Regularization 
(MPR) method proposed in |MG12| . and we determine pri- 
ors which follow the smooth phase-space structures traced 
by the weight distribution, as the latter adapts to match the 
observational data. This method facilitates recovering both 
a smoother and more accurate mass distribution function 
from noisy data, smoothing over local fluctuations without 
erasing global phase-space gradients. 



Maximizing the function ([2]) translates into a prescrip- 
tion, the so-called "force-of-change" equation, for correcting 
the weights of the particles while these are evolved in the 
total gravitational potential. In generic form this is 



dwi{t) 
dt 



-£Wi{t) 



A; (t) 



5 ; d k?\v. 

id 2=1 '■> 

dS . 9 ^, r 

3=1 



(4) 



The fi rst term on the right is the x 2 derivative term as in 
IDL07I . summed over the observational constraints for the dif- 
ferent types of data observables, id: Ai m , slits, and slitlets; 
Jid are the respective number of constraints. Aj and of 
axe the temporally-smoothed relative deviation and error in 
observable (id, j), and K 1 ? the corresponding kernel eval- 
uated at phase-space position Zi(t) of particle i. The last 
te rm ha s been written as sum over the log likelihoods as 
m iDLOSl . The notation is explained further in these papers. 
We have included here additional separate weight factors £id 
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for the different observables, which we use to optimize the 
influence of the da ta on the evolution of the model (see also 
lLong fe Macll201Ch . 



3.2 The gravitational potential 

The particle modelling technique allows for using both a 
fixed potential, known a priori, and a time- varying potential, 
self-consistently computed from the particle distribution. In 
our dynamical models, we assume that the total gravita- 
tional potential is generated by the luminous and dark mat- 
ter distributions, i.e. 

0= (5) 

Following ISellwoodl i|2003ft and IDL071 0* is frequently 
computed from the iV-particle model for the light distri- 
bution via a spherical harmonic decomposition, assuming a 
constant mass-to-light ratio T. The value of T is not a fixed 
parameter, but rather it is determined during the NMAGIC 
run, simu ltaneo usly with the modelling of the observational 
data fsee lDLOSl) . 

Instead, the dark matter halo potential <j)D is 
parametrized, and has the logarithmic form 



(R,z) = ^-hx(r 2 +R 2 



(6) 



l|Binnev fc Trem aine 200^), where vo and ro are a charac- 
teristic (constant) circular velocity and scale-length. This 
mass model ha s been widely and successfully used to fit 
galaxies (e.g. | Fall fc Efstathioul Il980l; iPersic et all [l99rj; 



iKronawitter et al 



20001 ; iThomas et al.ll2007bl . iDLOSh . 
Our dynamical models will explore a range of circular 
velocity curves, whose behaviour at large radii varies be- 
tween the near-Keplerian decline (when stars dominate the 
total potential), and a nearly flat (quasi-isothermal) shape 
obtained for massive dark halos, as shown in Fig. [15] 



3.3 The initial particle model 

We set up initial models of N = 750000 particles extending 
to 30.Re, i.e. ~ 110 kpc at the distance of NGC 4494. 

The density of these initial particle models is given 
by a spherical deprojection of the circularly- averaged sur- 
face brightness profile. We follow the method described in 
iGerhardl (|l99ll ) to obtain an isotropic stellar velocity distri- 
bution in the total gravitational potential generated by the 
stars plus the different dark matter halos described above. 

The particles' coordinates and velocities are drawn 
according to t he complete distribution func tion applying 
the method of iDebattista fc Sellwoodl l|2000r i. and particle 
weights are set equal to 1/N. 

Finally, following iKalnais! i| 19771 ) and |PL08l. we switch 
the sign of the velocity of a fraction of the retrograde parti- 
cles, with probability 



p(L z ) =p - 



L - 



(7) 



L'i + Lq 

(where po = 0.3, Lq = 0.02), to introduce some angular mo- 
mentum about the z-axis, while maintaining a smooth DF 
in equilibrium. This expedient makes it easier to reproduce 
the rotation velocity seen in the kinematic data. 



E 




Figure 6. The shaded area shows the range of circular velocity 
curves corresponding to the total gravitational potentials used 
in the dynamical models. Stars represent the self-consistent case 
with constant M/L; solid lines correspond to models embedded 
in various spherical dark matter halos. Black lines correspond to 
dark halos with ro = 4ii c and vo = [70, ...,270] kms -1 (from 
bottom to top). Solid and dashed red lines correspond to vq = 
150kms _1 , and rp = 3R C and 5R C , respectively. 



3.4 Photometric and kinematic observables for 
the modelling 

We now explain how the different photometric and kine- 
matic observational data are processed in order to be used 
as constraints for the NMAGIC models. 

For the photometric observables, the 3D luminosity den- 
sity profile obtained from the deprojection of the surface- 
brightness is expanded in spherical harmonic functions, and 
the expansio n coef ficients Ai m are used as luminosity con- 
straints fsee |PL07h . These Ai m are computed on a grid of 
50 quasi-logarithmically spaced radial shells between r m i n = 
0.01" and r max = 20R C , and all terms that would allow non- 
axisymmetry are set to zero. Poissonian errors are used for 
the mass in shells, whereas the errors for the higher-order 
mass moments are determined via Monte C arlo sim ulations 
of the density field of the target galaxy (see lDL07T ). 

For the long-slits and slitlets data, we use luminosity- 
weighted kinematic observables. Hence, we add to the tar- 
get observables the luminosity in each slit cell or slitlet, 
obtained by integrating the surface brightness distribution 
with a Monte Carlo algorithm. For the slits, we assume that 
the slitwidth equals 5"; for the slitlets, we set the slitwidth 
equal to the larger value between 1" and the diameter of a 
circle containing at least 250 particles of the initial particle 
model, so as to limit particle noise when computing model 
observables. These slit cell luminosities are independently 
fitted, with errors set to 1% of the luminosity in each cell. 

PNe data are modelled by maximizing the likelihood of 
the sa mple of discrete velocities and positions, as detailed in 
iDLOSl For computing the likelihood, the particles and PNe 
are binned in elliptical segments, assuming an average pro- 
jected ellipticity of 0.2. We consider 3 radial and 4 equally- 
spaced angular bins, with the first angular bin centered on 
the major axis, as shown in Fig. [7] Each segment contains 
at least 30 PNe. 
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Figure 7. Segments in which the line-of-sight velocity distribu- 
tion of particles is computed, for the likelihood method used to 
fit PNe with NMAGIC. 

The number of independent data constraints for the dif- 
ferent observables are JA lm = 450, Jsi = 2390, JLets = 346, 
and Jpn = 257. Since the errors are smallest for the Ai m , 
while the Poisson errors for the PNe are large, the Ai m 
terms contribute most to the force-of-change equation, while 
changes to the particle weights due the PNe data are small. 
To ensure that the halo is appropriately modelled, we in- 
crease the weight of the slitlets in the force-of-change by a 
factor of 4, and that of the PNe by a factor of 20. In practise 
we use the 4-folded data set for the slitlets ( 32.2.20 and the 
2-folded data set for the PNe ( H2.2.3[l . so the actual weight 
factors in equation Q are ^A lm = £si = ^Lcts = 1, and 
fp N = 10. 

3.5 Modelling procedure 

Starting from the initial particle model described in Sec- 
tion !3.3l the weights of all particles are evolved by NMAGIC 
until the A^-body system matches the target. During the 
whole evolution, the potential of the dark halo is kept fixed, 
but the stellar potential is frequently computed from the 
light distribution of the particle system. Particles are inte- 
grated in the total gravitational potential using a leapfrog 
scheme with adaptive time step. 

After a relaxation phase of 1000 steps in which the ini- 
tial particle model is advanced without any weight correc- 
tion, particle weights are updated according to the force-of- 
change obtained maximizing equation ©, for ~ 10 5 correc- 
tion time steps. Finally, the particles are freely evolved for 
another 10 4 steps without any further weight correction, to 
ensure that the final particle model is well phase-mixed. For 
reference, 10 4 correction time steps correspond to ~ 300 cir- 
cular rotation periods at the target R a . We define the model 
to have converged if % 2 averaged over 50 steps stays almost 
constant in the last 10 4 steps, with fluctuations which are 
typically of order 2%. 

During the weight adaptation, the model is regularized 
by the MPR method proposed in |MG12| . New p riors are de- 
termined in phase-space as described in IMG12I . by binning 
particles according to their orbital integrals. In particular, 
we grid particles according to energy (E), total circularity 



x, and angular momentum with respect to the rotation axis 
L z , in a grid of «b = 20, n x = 4, n z = 2 bins, which at 
the same time resolves the relevant phase-space structures 
and ensures enough particles in each grid cell. The average 
weight contained in each cell is computed, and then a thin- 
plate smoothing spline is fitted to the grid of priors, before 
assigning them to the particles. The individual particle pri- 
ors computed in this way are not kept constant in time but 
rather they are frequently updated while the particle weights 
are adapted to match the target observables. 



4 PARAMETER ESTIMATION WITH 

NMAGIC: HOW WELL CAN THE DARK 
MATTER HALO BE RECOVERED? 

In this work, we will be fitting the observational data of 
NGC 4494 with a sequence of NMAGIC models obtained 
for different dark matter halos and for different inclinations. 
In doing this, we wish to quantify the uncertainties in the 
best-fit parameters from a statistical distribution appropri- 
ate for our data and modelling method. In this Section, we 
describe a simple Monte Carlo method to estimate such a 
distribution. 

In the following, the best-fit parameters are defined as 
those for which a minimum is achieved in the merit function 

G =\x- A (8) 

which measures the agreement between the independent ob- 
servational data constraints and the final model in terms 
of the total y 2 from Ai m , slits, and slitlets that the x 2 M2M 
method tries to minimize, and of the log likelihood C for the 
PNe that the method tries to maximize. In order to carry 
out a robust parameter estimation, we also need to quantify 
the uncertainties in these best-fit parameters. 

Boundaries of the confidence regions for the estimated 
parameters will correspond to contours of constant AG rel- 
ative to the best-fitting model in the sampled parameter 
space. In common-practice dynamical mass analysis (see In- 
troduction), uncertainties in the model parameters are typ- 
ically estimated by considering the parameter change nec- 
essary to increa se \ 2 b y 1 (or 2.3), and to decrease the log 
likelihood £ (see lDLOSh by a factor of 0.5 (or 1.2) for 1 (or 2) 
fitted potential parameters (|Press et al.l 1992V This assumes 
Gaussian data, a model linear in its parameters, and that 
the number of fitted weight parameters N minus smooth- 
ing constraints N 3 is less than the number of data points J 
so that the number of degrees-of-freedom A*d f is positive. 
In fact, the dependence on the potential parameters is non- 
linear, and in our models N 3> J so that Adof may be less 
than zero. In addition, we find empirically for our data set 
and modelling set-up that A% 2 values are 3> 1. Therefore 
in this Section we compute confidence regions for parameter 
estimation via Monte Carlo experiments. 

In Section 4.1, we build a mock target galaxy closely 
similar to the elliptical galaxy NGC 4494, which is embedded 
in a dark matter halo with scale radius and circular velocity 
corresponding to the best-fitting model determined by N09 
for NGC 4494. We generate photometric and kinematic ob- 
servables for this mock galaxy with the same observational 
errors as for NGC 4494. In Section 4.2, we use the pseudo- 
data and known intrinsic parameters of this NGC 4494-like 
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target galaxy to calibrate the regularization parameter ap- 
propriate for the modelling of NGC 4494. In Section 4.3, we 
determine how well the circular velocity and scale radius of 
the known dark halo, the total enclosed mass, and the mass- 
to- light ratio of the target galaxy can be recovered, in a more 
rigorous way than in previous NMAGIC modelling (|DL08l ; 
iDLOStlPas et alj|201ll ). Finally, in Section 4.4 we use Monte 
Carlo experiments with such NGC 4494-like target galaxies 
to determine the appropriate AG to be used for confidence 
boundaries. This will then allow us to quantify uncertainties 
in the dark halo mass distributions of the real NGC 4494 in 
Section [5] 
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4.1 An NGC 4494-like galaxy and its observables 

Our initial NGC 4494-like target galaxy has the luminosity 
distribution obtained by deprojecting the surface brightness 
of NGC 4494 for i = 90° , and it is embedded in a logarith- 
mic dark matter halo (equation JBJ, with ro/R e = 4 and 
wo = 150kms _1 ), as in the best-fitting Jeans model of |N09| 
for NGC 4494. Similar to the orbital anisotropy predicted 
by so me merger models in the current cosmological scenario 
(e.g. iDekel et al.l 120051 ). the velocity distribution of this 
NGC 4494-like galaxy is isotropic in the center and increas- 
ingly radially anisotropic at larger radii. As for NGC 4494, 
we observe the target galaxy from a distance of 15.8 Mpc, 
and the projected effective radius R e ss 49". We set its stel- 
lar mass-to-light ratio to 3.8. 

Following the procedure outlined in Section |3]3j we gen- 
erate a spherical particle model realization for this target 
with 750000 particles. To implement the orbital anisotropy 
of the target galaxy, we a dopt a m o del w ith specified cir- 
cularity function following iGerhardl (|l99lf ) . We also add a 
certain amount of rotation to the particle model (see Sec- 
tion l3,3[l so as to mimic the real NGC 4494. With this spher- 
ical model as starting point, we then use NMAGIC to evolve 
it to the edge-on luminosity distribution of NGC 4494. The 
resulting axisymmetric model galaxy is then used to com- 
pute mock observables. 

As luminosity observables we use the coefficients Aim 
of the spherical harmonic expansion of the 3D luminosity 
density, computed as described in Section 13.41 The kine- 
matic observables are u, cr, /13, /u, projected onto both the 
long-slit set-up and the slitlets set-up used for the mod- 
elling of NGC 4494. These observables can be readily deter- 
mined using NMAGIC with weight correction turned off, in- 
tegrating the particle model fo r the ta rget galaxy, and time- 
smoothing the observables (see lMGlih . Observational errors 
from NGC 4494 are adopted, and Gaussian random variates 
with la equal to these errors are added to the kinematic 
observables computed in this way. 

Finally, we generate a mock sample of PNe velocities 
similar to that described in Section 12.2.31 We again use 
NMAGIC to integrate the particle model which imitates 
NGC 4494, and in parallel compute the time-averaged veloc- 
ity and velocity dispersion of the particles binned in radial 
and angular segments on the sky. We adopt the same 3 ra- 
dial bins and 4 equally-spaced angular bins as in Section ^. 41 
Then, we consider the catalogue of PNe used for the mod- 
elling of NGC 4494, and assign a new velocity to every PN 
in the catalogue, according to a Gaussian distribution with 
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Figure 8. Quality parameters for the final NMAGIC particle 
models as a function of the regularization parameter fi, from un- 
smoothed models (small values of fi) to oversmoothed models 
(high values of /i). Top panel: rms deviation (%) of the inter- 
nal velocity moments of the final particle model from the tar- 
get, computed over the radial region constrained by the data. 
Middle panel: merit function G in equation (BJ. Bottom panel: 
X 2 /J of the particle model fit to the J = 3186 photometric and 
absorption-line kinematic data points. 



velocity and velocity dispersion of the spatial segment the 
PN would belong to, given its position on the sky. 



4.2 Calibrating regularization 

The regularization parameter to be used when construct- 
ing dynamical models is case-dependent, and is influ- 
enced by several factors, most notably the observational 
data to be modelled (error bars, scatter, spatial cover- 
age) and the phase-spac e structure of the target galaxy 
(see iGerhard et al.l 1 19981; Cretton fc van den Boschl 1 19991 ; 
iThomas et al.ll2005f . lDL0a ; lDL09l ; lMG12l V 

To determine the optimal amount of smoothing given 
the data, we model the o bservat ional data of the NGC 4494- 
like galaxy using MPR (|MG12l ). varying the regularization 
parameter fi which controls the balance between regulariza- 
tion and goodness-of-fit in equation @. 

The results of our experiments are summarized in Fig.[8l 
which shows, for increasing values of fx, the reduced x 2 , the 
merit function G in equation (J5J, and the rms difference 
(Akin) between the internal velocity momen ts of the target 
and of the final NMAGIC model fsee lMG12l ). 

The minimum in the A^m plot determines the value of 
fj for which the model best recovers the internal moments of 
the input model. This occurs at fj, ~ 10 4 . In Fig. [8] we also 
see that for values of /x ^ 10 5 a smooth model fitting the 
data points within errors can no longer be found. 

Based on the results of these experiments, in the fol- 
lowing we will adopt a value of fi = 10 4 to regularize our 
NMAGIC models. 
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Figure 9. Recovering the halo parameters of NGC 4494-like 
galaxies. Each circle represents a final NMAGIC model. Top 
panel: the true dark halo (green triangle) has scale radius ro = 
4i? e and circular velocity no = 150kms — 1 . Bottom panel: this 
more massive dark halo (green triangle) has scale radius ro = li? e 
and circular velocity vq = 170kms — The colour scheme reflects 
the magnitude of the AG of each model relative to the overall 
best-fitting model, defined as that with the smallest G (red cross). 
The black lines correspond to 70% and 90% confidence contours, 
as derived from the simulations of Section I4.4I 



4.3 Recovery of input parameters for the dark 
halo of an NGC 4494-like galaxy 

We now construct NMAGIC models that fit the data of the 
NGC 4494-like target galaxy for a range of assumed pa- 
rameters of the dark halo potential in equation (J6]). The key 
question that we are interested in is: how well can NMAGIC 
recover the true parameters of the dark halo given the ob- 
servational data at hand? 

We follow the modelling procedure outlined in Sec- 
tion I3.51 starting from an isotropic initial particle model 
computed in each dark halo, as explained in Section I3.3I 
We use all available photometric and kinematic con- 
straints, including PNe through the likelihood technique. 
Our NMAGIC models are regularized with MPR, adopting 
the optimal value of jj, = 10 4 determined above. 

We sample the halo parameters (ro,t>o) on a grid of ~ 
200 models. The results of these experiments are shown in 
the top panel of Fig. [9] in terms of the AG of each NMAGIC 
model relative to the best-fitting model, which corresponds 
to the overall smallest G. The shape of the AG contours is 
regular, and there is an extended region for which the values 
of G are similar. Many models for different parameters of 
the dark halo provide similarly good fits to the data, all 
with y 2 / J ^ 1 (the number of data points entering \ 2 i s 
J = 3186). 

The best-fitting model is found for values of the dark 
halo parameters ro = 5.5i? e and vq = 180kms _1 , and it has 
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Figure 10. Total enclosed mass within 2/? c (stars plus dark mat- 
ter) for different values of the halo parameters ro and vo • Contour 
values are in units of 10 10 -Mq. 



X 2 /J = 0.22, log likelihood C = -770.5, and G = 1125.6. 
The scale radius of this best-fitting model is larger than that 
of the true model, which achieves \ 2 1 J = 0.24, log likelihood 
C — —769.9, and G = 1148.7. For comparison, the no-dark 
halo model, which provides the poorest fit to the available 
observational data, has a final xyJ = 0.37, log likelihood 
£ = -782.2, and G = 1376.6. 

Surprisingly, the values of AG are large (3> 1). For in- 
stance, the model obtained for the true parameters of the 
(known) dark matter halo achieves a value of G which is 
AG ~ 20 (Ax 2 ~ 60) away from the best-fitting model. 
Applying the classical \ 2 statistics for two degrees of free- 
dom, the true parameters of the dark matter halo would be 
discarded at the 99.99% level. We investigate this issue in 
Section ECU 

Note that the contours of AG re main open to the top- 
right edge of the plot (see also e. g. | Gerhard et all 1 19981 ; 
iThomas et"aill2005l ; iMurphv et al.ll201lF 7 showing that it is 
not possible to put robust constraints on both halo param- 
eters simultaneously in this case. Similar values of G, and 
hence similarly good models are achieved for halo models 
located along a large band extending from low ro and low 
vo to high ro and high vo- 

This diagonal band shrinks when modeling target galax- 
ies embedded in more massive halos, as the models that do 
not contain enough mass are then ruled out. This is shown 
in the bottom panel of Fig. [9j which is obtained by repeat- 
ing the analysis for an NGC 4494-like galaxy embedded in 
a more massive dark matter halo (whose parameters are set 
to ro = lR e and vo = ^Okms^ 1 ). 

Models inside the diagonal band share a similar total 
enclosed mass within the radial region constrained by the 
observational data, as shown in Fig. [TD] by means of con- 
tours of total mass within 2R C in the same (ro,«o) parame- 
ter space. Evidently, the total enclosed mass inside a certain 
radius is what the dynamical models constrain best with the 
data at hand. Indeed, when plotting the values of AG of the 
final NMAGIC models against enclosed mass (see Fig. Ill|l 
the shape of a parabola is evident, despite the scatter in the 
values of AG. 
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Figure 11. AG of the final NMAGIC models fitting the data of 
the NGC 4494-like galaxy, as a function of the total mass inside 
2R B divided by the corresponding quantity for the best-fitting 
model. The red cross shows the best-fitting model, the green tri- 
angle the model with the true dark halo parameters (ro = 4i? c , 
vq = 150kms -1 ). The black lines correspond to 70% and 90% 
confidence regions derived from the simulations of Section 14.41 



Figure 13. Recovery of the dark matter fraction /dm = A^dm/-^ 
as a function of radius for the range of explored NMAGIC models 
of the NGC 4494-like data (orange shaded region). The green 
line represents the known target, and the red line the best- fitting 
model. The yellow shaded region shows the 70% confidence region 
(see Section l4.4l l. 



200 - ' ' 



150 - 




3.2 3.4 3.6 3.8 4.0 

M/L 



Figure 12. As Fig. [TT] but AG of the final NMAGIC models 
as a function of stellar mass-to-light. The galaxy has (known) 
M/L = 3.8. 



Figs. [P2lfl"3l show how well the NMAGIC models re- 
cover the mass-to- light ratio and dark matter fraction of the 
NGC 4494-like galaxy from the mock data. These are both 
important quantities which we are interested in measuring 
with accuracy for real galaxies. 

At this stage, two obvious questions are: what is the typ- 
ical AG difference between the best-fitting NMAGIC model 
and the model with the true halo parameters? And related 
to this, within what errors can we trust the dark halo pa- 
rameters of the best-fitting model? The following analysis is 
designed to answer these questions. 
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Figure 14. Normalized cumulative distribution of AG (left) and 
Ax 2 (right) between the true dark halo and the best fitting model, 
derived from NMAGIC models of the NGC 4494-like galaxy em- 
bedded in 10 different dark halos. The best fitting model is defined 
as that with the minimum G (x 2 )- The red lines correspond to 
70% and 90% confidence levels. 



4.4 Parameter estimation for the dark matter 
halo: confidence levels 

We use Monte Carlo simulations to estimate the values of 
AG which represent a specific confidence level, given the ob- 
servational data of NGC 4494 and our NMAGIC modelling 
technique. 

In practice, given constraints on computer time, we con- 
struct a sequence of NGC 4494-like galaxies embedded in 10 
different dark halos, and model each one of them with a 
range of dark halo parameters. Then, for each target galaxy 
we compute the AG between the best-fitting model and the 
model obtained for the true parameters. In this way, we can 
estimate the values of AG within which the true parameters 
are found in 70% or 90% of our experiments. 

The main results from these experiments are summa- 
rized in Fig. 1141 which shows the cumulative distributions 
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of AG and A\ 2 between the true dark halo model and the 
best fitting model, for the 10 different experiments. It is seen 
that the NMAGIC models recover the true potential of the 
dark matter halo within AG < 26 about 70% of the time, 
and within AG < 32 about 90% of the time. Corresponding 
values when the log likelihood is not included are: A\ 2 ^ 59 
and Ax 2 < 70 for 70% and 90% confidence regions, respec- 
tively. 

The magnitude of these differences is quite surpris- 
ing. Therefore we perform several additional tests to un- 
derstand them better. First, we verify that the measured 
differences AG are actually significant with respect to fluc- 
tuations caused by modelling noise or measurements uncer- 
tainties. Indeed, it is natural to speculate that the numerical 
noise in the procedure of the weight-adjustment may cause 
fluctuations of G over time. However, these fluctuations are 
comparatively small; on average they are ~ 8 in G (~ 17 in 
X 2 ), i.e. a factor of three smaller than the typical values of 
AG (and A x 2 ). 

Then, in order to quantify the relevance of the obser- 
vational measurement uncertainties, we generate a Monte 
Carlo chain of one hundred realizations of the data of the 
NGC 4494-like mock galaxy of Section ft.ll Every realization 
is constructed by drawing random Gaussian-distributed val- 
ues for all kinematic data points, such that the mean is as 
predicted for the mock galaxy, and the variance corresponds 



2500 



to the observationa l erro r s (see e.g. van der Marel et alj 



1 19981 ; ICretton et ail |2000| ; I Thomas et afll2005h . Then, we 



model each realization with NMAGIC, assuming the true 
gravitational potential of the target galaxy. Based on mod- 
eling 100 realizations of the observational data, we conclude 
that random variations in the data according to the observa- 
tional errors correspond to fluctuations of ~ 16 in the merit 
function G (and ~ 32 in \ 2 ) within la. Instead, the exper- 
iments above showed that if the potential is not known in 
advance, NMAGIC models match the true potential of the 
dark matter halo within AG < 26 at la level. Hence, the 
A\ 2 differe nce derived from only fitting data sets of the true 
model (e.g. iThomas et al.ll2005l ) appears to be an underes- 
timate. 

As a further check, we verify whether the minimum 
value of G and/or the magnitude of AG are influenced by 
the degree of regularization employed in the NMAGIC fits. 
In the framework of Schwarzschild modelling, IVerolme et alj 
(2002) showed that smoothing constraints do not affect 
the shape of x 2 -contours significantl y, and leave th e over - 
all best-fitting parameters unchanged. IThomas et all (|2005l ) 
constructed Schwarzschild models for different values of 
the smoothing parameter, and computed the A^ 2 due to 
-/Vmock = 60 different realizations of the data sets as a 
function of the smoothing parameter. They found that the 
magnitude of A\ 2 increases for increasing regularization 
(see their Fig. 6). We repeat their experiment modelling 
Nmock = 100 realizations of the data of the NGC 4494-like 
mock galaxy for different values of the regularization pa- 
rameter /i, and show our results in Fig. 1151 The error bars 
represent the la-variation AG as a function of /i (note that 
they are multiplied by a factor 5), and show no indication 
that such AG is influenced by the regularization, up until 
and including the optimal regularization parameter used in 
our modeling (fj, = 10 4 ). An increase of the fluctuations of G 
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Figure 15. The merit function G = x 2 /2 — C for increasing values 
of the regularization parameter /i. Dots represent the mean over 
100 model fits to Monte Carlo realizations of the kinematic data of 
the NGC 4494-like galaxy; error bars represent the rms deviation 
from the mean multiplied by a factor 5. 



due to regularization is only seen for oversmoothed models 
(jj, = 10 5 , see Section l4~2j). 

The Monte Carlo experiments also reveal some slight 
biases which are intrinsic to our diagnostics of the best- 
fitting model. In particular, we find that the minimum \ 2 
values are achieved on average for slightly more massive ha- 
los than that of the target galaxy, as previously noted in 
iGerhard et"afl \ 19981 ) . while the maximum C values are typ- 
ically achieved for halos which are slightly more diffuse than 
the true (known) halo. The combination of \ 2 an( i C in the 
merit function G in equation |(8j is less biased, since it strikes 
a compromise between these two opposite trends. 



4.5 Summary 

To summarize, we investigated how well the dark halo pa- 
rameters can be constrained from data which have the spa- 
tial coverage and error bars of the current observational 
data for NGC 4494. These data do not suffice to determine 
uniquely both the scale radius and the circular velocity of 
the halo specified in equation ([6]), and different combina- 
tions of (ro,vo,T) that yield similar values of the total en- 
closed mass are allowed. The enclosed mass within 2R e can 
be determined to within 10%, and the dark matter fraction 
/dm(3-R c ) to within ±0.1. For the family of logarithmic dark 
matter halos, the total circular velocity can be determined 
to ±20kms _1 . 

Via Monte Carlo experiments, we computed confi- 
dence levels for parameter estimation. We conclude that, 
for NMAGIC models of the observational data at hand, the 
usual la (70%) level corresponds to a value of AG ~ 26, and 
of Ax 2 ~ 59 (and AG ~ 32, A x 2 ~ 70 at 90% confidence 
level). Using these statistical results, we can now derive the 
uncertainties on the determination in the halo parameters 
for NGC 4494. 



© 0000 RAS, MNRAS 000, 000-000 



Dark matter halo of NGC UH 13 



5 DYNAMICAL MODELS OF NGC 4494 

We now use NMAGIC to construct axisymmetric dynam- 
ical models for NGC 4494 fitting all the photometric and 
kinematic data described in Section [2] Three different in- 
clinations for the stellar distribution are considered, i = 
90°, 70°, 45°, for which we carried out the deprojection of 
the photometric data (see Section l2.ip . The latter value is 
close to the minimum i nclination allowed by the observed 
flattening of NGC 4494 (|N09| ; IF1 ll ). For each inclination, we 
explore a sequence of gravitational potentials, including the 
self-consistent case with constant M/L and various quasi- 
isothermal dark matter halos. 

The main results of our suite of dynamical models are 
presented in Fig. 1161 In analogy with the analysis of the 
previous Section, we plot the AG of each NMAGIC model 
relative to the best-fitting model, defined as that with mini- 
mum G, for each inclination i — 90°, 70°, 45°. For reference, 
the best-fitting model for i = 90° has \ 2 1 J = 0.41, log like- 
lihood C = —793.3, and G = 1449.2. In the same figure, 
we also overplot the 70% and 90% confidence levels deter- 
mined from the Monte Carlo experiments of Section f4. 41 In 
the right panel, we show the final mass-to-light ratios of our 
dynamical models. 

The shape of the AG-contours is regular, and the range 
of dark matter halos consistent with the data within the 
confidence levels has circular velocity in [160-240] kms -1 
and scale radius ~ l-2i? c . The stellar mass-to- light ratio is 
Ty ~ 3. The main characteristics of the mass distribution 
of the explored models, and of the preferred models, are 
plotted in Fig. 1171 The total circular velocity curve of the 
good models is approximately flat ("isothermal") outside 
1R C , with v c {3R c ) ~ 220kms~ 1 . Less massive halos, as well 
as models with constant M/L, are not consistent with the 
data. 

It is apparent from Fig. [16] that these main results, and 
the topology of the contours of G, are independent of the 
adopted inclination of the stellar distribution of NGC 4494. 
In general, we find that edge-on models (i = 90°) are pre- 
ferred, in the sense of a lower value of G, a lower value of 
X 2 , and a higher value of C, respectively. A discussion of this 
issue is deferred to Section[5]3] whereas we now describe the 
model fits to the data, and the orbital structure of the dy- 
namical models. 



5.1 Model fits to the observational data 

For a wide range of dark halos, NMAGIC finds very good 
fits to the observational data in terms of \ 2 values. In par- 
ticular, for all the assumed values of the inclinations, the 
models embedded in the dark matter halos that that we 
tried converged to \ 2 1 J values per data point less than 1 
(except for extreme values of the parameters, corresponding 
to the bottom-right corner of Fig. I16|l . The models with con- 
stant M/L achieve the poorest fit, with a value of reduced 



0.60 (90°), x /J = 0.94 (70°), and x /J = 1-15 



x 2 /J 

(45°). 

The fits to the photometric constraints are generally 
excellent, and they are visually indistinguishable for most 
potentials, and for the different inclinations. Models com- 
patible with the data at the 90% level achieve these fits 
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Figure 16. Results of the NMAGIC dynamical models of 
NGC 4494 for the range of explored halo parameters (ro,«o). From 
top to bottom: i = 90° , 70° , 45° . The grey contours correspond to 
70% and 90% confidence levels, as determined in Section l4.4l Left 
column: AG of NMAGIC models relative to the best-fitting model 
(separate colour bar for each inclination). The bottom-right cor- 
ner is a region in which no good models for the data could be 
found. Right column: V-band mass-to-light Ty of the final mod- 
els. 



i — 90° , 70° , 45° , respectively. Fig. [18] shows the very good 
fit to the first moment of the Ai m , i.e. the differential stellar 
mass distribution, for the preferred inclination of 90° . Three 
fiducial halo models are shown: the best-fitting model, and 
two models which fit the data with x 2 < li but ue outside 
the 90% confidence band, having a more massive and a more 
diffuse dark halo than the best-fitting models. 

The fits to the long-slit kinematics are also gener- 
ally good. Models within the 90% confidence range achieve 
a maximum value of X 2 1 Js\ = 0.36,0.41,0.35, for i — 
90°, 70°, 45°, respectively. For the fiducial models, these fits 
are shown in Fig. 1191 where it is apparent that the cen- 
tral feature of the kinematically decoupled core in v and in 
I13 is reproduced well by our particle models, and that the 
models have some difficulty in matching the detailed long- 
slit kinematics along the outer major axis. In particular, 
the best-fitting models have slightly lower values of v and 
slightly higher values of a. To a great extent, these system- 
atic deviations can be traced to the compromise that the 
models must find (see below) between the long-slit and the 
slitlets kinematic data, which extend to larger radii. The 
fiducial models show that more massive/diffuse halos result 
in higher/lower values of a and hi along both axes. Finally, 
the shaded regions in Fig. \T§\ show the range of models al- 
lowed by the different inclinations (for the more massive and 
the more diffuse fiducial models) , and by the 90% confidence 
region (for the best-fitting model). 

The fits to the newer slitlets kinematics are achieved 
with maximum xV^Lcts = 1., 1.18, 1.27 within the 90% con- 
fidence level, for the three inclinations of 90° , 70° , 45° , re- 
spectively. Fig. [2D] shows the fits of the preferred edge-on 
particle models to the slitlets kinematics for the fiducial po- 
tentials. The trends with halo mass described before are still 
clearly visible. Moreover, it can be seen that the best-fitting 
edge-on model, corresponding to the black line, has to find 
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Figure 17. Top /ine panels: as a function of radius, the total 
mass, the dark matter fraction /dm = J^dm/M within R, and 
the circular velocity curves corresponding to dark matter, stars, 
and total potential. The orange shaded region shows all edge-on 
NMAGIC models investigated, the yellow shaded region shows 
the 70% confidence region (see Section |4.4|| . Bottom panel: the 
contributions of stars (red) and dark matter (grey) to the total 
circular velocity curve (violet) for the edge-on models in the 70% 
confidence band. 




Figure 18. NMAGIC fits to the differential stellar mass distri- 
bution (Aoo coefficient) for an inclination of 90°. The black line 
is for the best-fitting model (ro = lR e , v$ = 190kms — 1 ), the red 
line is for a model residing in a more diffuse halo (ro = 5R C , 
vq = 90kms -1 ), and the green line represents a model obtained 
for a massive halo (ro = lRa,VQ = 250kms _1 ). The three model 
curves are nearly identical even though the last two models are 
outside the 90% confidence region. 
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Figure 19. NMAGIC fits to the long-slit kinematic data (black 
dots) of NGC 4494 along the major (top) and minor (bottom) 
axis. The black shaded region shows the best-fitting edge-on mod- 
els within the 90% confidence range. The red and green shaded 
regions show the spread of two fiducial models for the explored in- 
clinations; red is for the model residing in a lighter halo (ro = 5R C , 
vo = 90kms -1 ), and green is for the model obtained in a massive 
halo (ro = lR e , Vq = 250kms ). Both models are well outside 
the confidence range. 



a compromise between the kinematic data in slit and slitlets 
in the outer region. This results in a lower magnitude for the 
velocity, and in a higher velocity dispersion, which explains 
the deviations from the long-slit kinematics in the outermost 
data points. 

Finally, Fig. [21] shows a comparison of the final 
NMAGIC models with the PN kinematic data, for all the 
inclinations and for the fiducial potentials. The velocity, ve- 
locity dispersion, and LOSVD are plotted in the angular 
segments used for the modelling (see Section 12.2.30 . The 
axisymmetric nature of our models is apparent in the re- 
flection symmetry of the kinematics in diagonally opposite 
segments. While the mean velocity profiles are well fit by all 
models shown, the low (high) dark matter models are sys- 
tematically low (high) in the PN dispersion plot, and show 
systematic deviations from the data histograms. The best- 
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Figure 21. NMAGIC fits to the PNe velocities for all the considered inclinations. Top two rows: mean velocity dispersion and velocity 
profiles in the angular segments centered on 0° (major axis), 90°, 180°, 270°, from left to right. Bottom three rows: LOSVD in the same 
angular segments at radii of 52", 110", and 300", going upwards. Colours of the shaded regions are as in Fig. 1191 

fitting models, instead, provide a good match to all the PN 
data. 



5.2 The internal kinematics of NGC 4494 

We now look at the intrinsic kinematics of the NGC 4494 
models, in order to learn about the orbital anisotropy of this 
galaxy. Fig. [22] shows the internal kinematics for the range 
of explored inclinations. 

As expected, in the radial region well constrained by the 
kinematic data the ratio of the radial to tangential velocity 
dispersions is larger for models embedde d in a massive halo 
than for low dark ma tter models (see e.g. iBinnev fc Mamonl 
1 19821 ; iGerhardll 19931 ). 

On average, the ratio of the radial to tangential veloc- 
ity dispersions along R in the best-fitting models is ~ 1.2, 
i.e. the underlying orbital distribution is slightly radially bi- 
ased. As discussed further in Section [6] the degree of radial 
anisotropy of NGC 4494 is milder than what was previously 
found for t he two other intermedi ate- lum inosity ellipticals 
NGC 4697 i|DL08l ) and NGC 3379 (|DL09h . 

For all the explored potentials, the two components of 
the tangential velocity dispersion are similar along the mi- 
nor axis, as required for an axisymmetric system, whereas 
the azimuthal dispersions are higher in the equatorial 
plane, suggesting together with the radial anisotropy that 
this (rotating) el liptical may be flattened by meridional 
anisotropy (e.g. iDehnen fc Gerhard! 1 19931 ; iThomas et al.l 
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Figure 20. NMAGIC fits to the slitlets kinematics near the major 
axis of NGC 4494 (black dots with error bars). For comparison, 
light-green dots represent the original major-axis slit data. The 
model points are averages over the same slit cells as the data, 
and are connected by straight line segments. Different colours 
represent the fiducial models: black is for the best-fitting model 
(ro = 1-Ro, Vo = igOkms" 1 ), red is for the model residing in 
a lighter halo (ro = 5i? e , Vo = 90kms~ 1 ), and green is for the 
model obtained in a massive halo (ro = IRe, vo = 250kms — 1 ). 
For simplicity, only the preferred edge-on models are shown. 



© 0000 RAS, MNRAS 000, 000-000 



16 L. Morganti et al. 




Figure 22. Internal kinematics of the final NMAGIC models 
for NGC 4494. The ratio of the azimuthal to meridional velocity 
dispersion (top) and of the radial to tangential velocity dispersion 
<Tt = \j ( a \ "b^ST^ (bottom) are plotted as a function of major 
axis R in the equatorial plane (left column), and of minor axis 
z in the meridional plane (right column). Colours of the shaded 
regions are as for Figs. [T9l 1211 The vertical dashed lines show the 
extent of the kinematic data (long-slit, slitlets, and PNe, going 
outwards) . 



2009a; iDas et all l2010bh . Variations due to the assumed 
inclination are always larger in the ratio between the az- 
imuthal and meridional velocity dispersions along R. 

Finally, at larger radii, and close to the outermost data 
points (marked by the vertical dashed lines), the orbital 
structure is nearl y isotropic along both R and z. As we have 
shown in IMG12I . in those external regions, which are un- 
constrained by the observational data, a bias towards the 
dynamical structure of the initial (isotropic, in this case) 
particle model cannot be avoided. 



5.3 The inclination of NGC 4494 

From our NMAGIC models of NGC 4494, we would con- 
clude that the observational data prefer an inclination of 
90° out of the three inclinations that we explored. Indeed, 
edge-on models provide better fits to the data than mod- 
els with lower inclinations, in terms of the merit function 
G, as well as \ 2 an d C In particular, we find a AG ~ 144 
(78) between the best-fitting edge-on model and the best- 
fitting model obtained for 45° (70°), which is highly signif- 
icant when compared to the typical differences AG among 
different potentials for a given inclination. Visually, it is hard 
to distinguish the fit to the observables for the 90° and the 
70° inclinations, whereas the more face-on 45° models are 
characterized by a systematically low velocity and velocity 
dispersion with respect to the observational data. 

Should this preference for edge-on inclinations really be 
believed? Up to now, we have assumed a spherical dark halo 
for all the considered inclinations, and only the stellar dis- 
tribution was axisymmetric and varied according to the in- 
clination. We can relax the assumption of sphericity of the 
dark matter halo, and investigate dark matter halos with 
similar flattening as for the stars. To this end, we consider 
the flattened logarithmic potential 
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Table 1. Results of NMAGIC fits to an NGC 4494-like galaxy 
whose stellar distribution has a known inclination iq, assuming 
an inclination im in the modelling. \ 2 values are normalized by 
the number of observables J = 3186. 



i>B(-R,2) = yln 



r'o + R 2 + 



(9) 



(|Binnev fc Tremaine|[200ci ). Then, we use the approximate 
relation 1 — ~ (1 — q p )/3 to derive the flattening of the 
potential q$ from that of the density distribution q p , which 
we infer from the apparent flattening q of the isophotes via 



2 • . 2.2. 

cos i + q p sin i. 



(10) 



For an average flattening of the isophotes q — 0.82, 
the intrinsic flattening of the density distribution q p — 
0.82,0.79,0.59 for i = 90°, 70°, 45°, respectively, and the 
corresponding q$ = 0.94, 0.93, 0.86. 

We have explored a few flattened dark halos with these 
parameters and the same spherically-averaged mass distri- 
bution as the best-fitting models for i = 70° and 45°. Al- 
though some of them provide slightly better fits to the obser- 
vational data, our main result does not change, and edge-on 
models are still preferred by far. 

The obvious next question is whether there could be 
some bias in the modelling techniq ue. This explanation was 



proposed bv lThomas et al 



(2007a), who found that all their 
best-fitting Schwarzschild models of A r -body merger rem- 
nants have i — 90°, and argued that edge-on models neces- 
sarily have a smaller \ 2 than face-on models, due to their 
greater freedom in the adjustment of prograde and retro- 
grade orbits to fit rotation and axisymmetric deviations from 
a Gaussian LOSVD. 

We used our NGC 4494-like galaxy (see Section I4.1jl 
to test whether particle models exhibit a bias analogous 
to that of orbit-based models. First, we constructed a se- 
ries of NGC 4494-like galaxies assuming different inclina- 
tions ic, = 90° , 70° , 45° for the stellar distribution. Then, we 
modelled the observational data of these mock galaxies with 
NMAGIC for different inclinations iyi- The results of these 
experiments are summarized in Table [H and con firm the 
findings and the argument of iThomas et all (|2007ah : on av- 
erage, the values of x 2 , and of the merit function G, increase 
for lower values of the inclination. Table [1] reveals that edge- 
on models generally provide better fits. This is true both in 
the case in which the model has the true known inclination 
of the target galaxy, and in the case in which the model as- 
sumes a wrong value of the inclination. For reasons that are 
currently unclear to us, the likelihood favours instead more 
face-on models. 

This test does not onl y confirm the existence of a bias 
(e.g., IThomas et "aLll2007al ). but also enables us to quantify 
its effect on the dynamical models of NGC 4494. The exper- 
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iments show that, when modelling a NGC 4494-like galaxy 
which has a true inclination ig = 45°, a AG ~ 95 between 
the edge-on model and the more inclined model could be as- 
sociated with the "edge-on bias". However, our NGC 4494 
runs display a significantly larger difference AG ~ 144 be- 
tween the edge-on and the im = 45° models. Thus, com- 
bining these results suggests that NGC 4494 is truly close 
to edge-on. If so, this might provide also a simple explana- 
tion for the slightly larger rotation velocities measured in 
the major axis slit compared to the surrounding slitlets (see 
Fig. I20[) . This could be the kinematic signature of a faint 
disk. 



6 DISCUSSION 

In this Section we discuss (?) the results on parameter esti- 
mation of dark matter halos with NMAGIC, (ii) our dynam- 
ical models for NGC 4494 in comparison with previous work, 
and (Hi) the dark matter distribution and orbital structure 
of NGC 4494 in the wider context of intermediate-luminosity 
elliptical galaxies with steeply falling, 'quasi-Keplerian' ve- 
locity dispersion profiles. 



6.1 Confidence limits for parameter estimation 
with NMAGIC 

In Section [4] we used Monte Carlo simulations of a 
NGC 4494-like mock galaxy in different dark matter haloes 
to assess the confidence limits for the halo parameters that 
provide statistically valid models for the real NGC 4494. 
The data set for the mock galaxy closely resembled that of 
NGC 4494. We estimate the differences in the merit func- 
tion AG corresponding to confidence levels enclosing a cer- 
tain probability of finding the true values o f the parameters 
(e.g. iPress et all 19921 ; iThomas et alj|2005h . Specifically, we 
constructed and modelled 10 NGC 4494-like galaxies in dif- 
ferent dark matter halos, and estimated the distribution of 
AG values between the best-fitting model and the model 
with the true halo parameters for our data and modelling 
set-up. 

These experiments showed that our NMAGIC dynam- 
ical models match the true potential of the dark matter 
halo within AG < 26 (Ax 2 < 59) about 70% of the time, 
and within AG < 32 (A X 2 < 70) about 90% of the time. 
These values also correspond approximately to the fluctua- 
tions found within the confidence boundaries. The numerical 
noise in G or x 2 caused by the adjustment of particle weights 
explains only a minor fraction of these AG or A^ 2 fluctua- 
tions. Fluctuations induced by varying the data within their 
error bars are much more significant compar ed to the mea- 
sured AG or Ax 2 values, as already noted bv lThomas et all 
|2005l ) in the context of Schwarzschild's method. Our ex- 
periments show, however, that some additional uncertainty 
must be related to the freedom associated with comparing 
different trial potentials. They also indicate that these re- 
sults are unaffected by the strength of regularization as long 
as the model is not oversmoothed. 

The relative differences AG (and A^ 2 ) found in our 
simulations are substantially larger than what is often as- 
sumed, and reveal that the A^ 2 = (1,2.3,..) criterion for 
v — (1, 2, ..) free model parameters used in many dynamical 



studies in the literature (see Introduction) is inappropriate 
for our NMAGIC particle models. 

The Ax 2 (^) approach in dynamical modelling is based 
on x 2 -statistics and assumes Gaussian errors, linear depen- 
dence on the model parameters, and that the num ber of 
degrees of freedom is positive (e.g. IPress et al.lll992j ). Then 
X 2 can be marginalized over all degrees of freedom other 
than the number of free parameters in the mass model. This 
implies the assumption that much of the freedom in the 
fitted weights (or orbits) is used by the regularization in 
constraining the distribution function from the data in the 
given trial potential. Because the number of particle weights 
is much larger than the number of data points, and given 
the indications that the simulations results are insensitive to 
the regularization parameter (Section 14. 4[) . this assumption 
could be incorrect. 

The A^ 2 values found in our experiments appear to be 
larger than in typical Schwarzs child applications. One excep- 
tion a ppears to be the work of lvan den Bosch fc van de Venl 
(2009) who determined confidence limits on the modelling 
parameters based on the ex pected standard deviation of \ 2 
itself, s(\ 2 ) = y2( J — M), where J is the number of ob- 
servational constraints and M = 4 the number of their free 
model parameters. However, the fact that the standard devi- 
ation of x is s (x ) does not invalidate the stand ard Ay (v) 
meth od when the underlying assumptions (e.g. IPress et al.l 
1992) are met. Thus, while the Ax 2 values found in our 
Monte Carlo simulations do have the same order of magni- 
tude as s(x 2 ), it is not clear whether this is significant until 
more simulations in different modelling contexts have been 
analyzed. 

The experience in dynamical modelling with made-to- 
measure particle models is still limited, and it is possi- 
ble that there are some aspects of these methods that in- 
fluence parameter estimation in different ways than, e.g., 
Schwarzschild models. However, given the wide-spread use 
of dynamical modelling for measuring, e.g., dark halo pa- 
rameters and black hole masses in galaxies, it is important 
to test the statistical premises of this work more thoroughly. 
Monte Carlo simulations like those we have performed may 
be the best way to tackle these issues, and determine the 
appropriate Ax 2 values for estimating confidence limits for 
a given observational and modelling set-up. 



6.2 Dynamical models for NGC 4494: comparison 
with the literature 

The elliptical galaxy NGC 4494 has previousl y attracted the 

van der Marell 



atten ti on of several dynamical s tu dies (e.g 



1991 



2001 



2011 



iKronawitter et al 



Romanowskv ct al 



|2000| ; iMagorrian fc Ballantvn. 



20031; iRodionov fc Athanassoul 



iNapolitano et al.ll2009l ; lLackner fc Ostrikedl2oToT ) 



Our NMAGIC models improve on the models explored 
so far in some important aspects. First, we considered as 
many observational data as currently avail able, i.e. at large 
radii we had availab le bo th PNe velocities (|N09f ) and galaxy 
spectra in slitlets (|Fllh . The newer slitlets data show a 
milder drop of the velocity disper sion than what it is sug- 
gested by the long-slit data IjFllI) . and this is likely to ac- 
count for the larger enclosed mass be yond ~ 2R a in our re- 
sults compared to the Jeans models of lN09l Second, because 
of the greater constraining power of the new data compared 
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to the PNe, we sampled the dark halo parameter space much 
more finely than before, and for different values of the in- 
clination (i — 45°, 70°, 90°). Thirdly, we performed a thor- 
ough analysis of the confidence levels at which these param- 
eters can be estimated with our models, given the data at 
hand. Finally, contrary to most previous studies that con- 
sider ed spherical models (but see lRodionov fc Athanassoulal 
2011), our NMAGIC models are axisymmetric, and we also 
explored the possibility of flattened dark halos. While the 
most robust results should eventually be derived using tri- 
axial models, it has been shown that relaxing the spherical 
assum ption hardly influences the recovered halo mass (see 
e.g. iDLOgf) . 

Although NGC 44 94 was among the three g alaxies de- 
scribed as "naked" in iRomanowskv et al.l (|2003l ) , and has 
an unusually l ow d ark matter fraction in the analysis of 
iDeason et all <|2012t ). our dynamical models show that a 
dark matter halo is require d . In th is respect, they agree with 
iRodionov fc Athanassoulal (|201ll ). although the small se- 
quence of explored models did not allow th em to put r obust 
constraints on the halo mass, and also with [N0!|. |N0^ inves- 
tigated a family of multi-component kurtosis-based Jeans 
models. Their favoured logarithmic dark halo had ro = 4_R e 
and vo = 150kms _1 . Our best-fitting halo models are more 
massive than what they obtained (see Section [5]), although 
probably consistent with their errors. These differences are 
probably mostly due to the additional slitlets constraints 
that we included, and the better use of the higher-order 
kinematic moments in the NMAGIC modelling. 

The dark matter fraction of our best-fitting NMAGIC 
models /dm(< 5R c ) = 0.6 ±0.1 is therefore higher than 
what was reported by |N09| (0.2-0.5). Interestingly, the new 
value agrees very well with the ste llar population predictions 
computed by IDeason" et al.1 (j2012r ) assuming a Chabrier ini- 
tial mass function (see their Fig. 7). Their spherical distri- 
bution function models of the PNe velocities, on the other 
hand, indicated a low dark matter fraction /dm(< 5R e ) = 
0.32 ±0.12. 

Alternative estimates of the mass of NGC 4494, e.g. 
from the X-ray emission of hot gas, would be highly desir- 
able. F rom the compact X-ray gas emission, iFukazawa et alj 
(2006) estimated Tb = 6.2 ± 1.9 inside lR e , consistent with 
the Jeans modelling of lN09l but higher than our NMAGIC 
results. However, the existence of X-r ay emitting gas around 
this galaxy has been qu estioned |0'Sullivan &: Ponmanl 
l2004lDiehl fc StatlerlkoOTl ), and the validity of hydrostati c 
equilibrium may also be dubious l|Ciotti &: Pellegrimll 19961 ). 

Additional kinematic constraints at large radii could be 
obtained from different tracers of the mass distribution, such 
as GCs. However, from the analysis of the spatial and kine- 
matic distributions of blue and red G Cs ar ound NGC 4494, 
which do not follow that of the stars (|Flll ). it is likely that 
PNe and GCs are in distinct dynam ical equilibria i n the 
same gravitational potential (see e.g. iDas et aT1l201lh . 

Our best-fitting NMAGIC models have a mass-to-light 
Tb ~ 3.71 ± 0.15, obta ined converting from T y using the 
de-reddened colour from lGoudfrooii et al.l ll 19941 ). This value 
of the stellar mass-to-light ratio is lower than that obtained 
from the Jeans models of[N09, and is easier to reconcile with 
independent measurements from stellar population models, 
whi ch gav e Tb — 4.3±0.7 for a Kroupa initial mass fu nction 
(see lN09f) . In this respect. lLackner fc Ostrikerl (|201Ch fitted 



galaxy formation model to the PNe velocities of NGC 4494, 
and found that the best-fit dissipational and dissipationless 
models give Tb = 2.97 and Tb — 3.96, respectively, so 
that a purely dissipation al formation scenario for N GC 4494 
seemed to be ruled out |Lackner Sz Ostrikerl |2010| ). 

Finally, the edge-on NMAGIC models that best repro- 
duce the observational data of NGC 4494 are mildly radially 
anisotropic, with (3 ~ 0.4 for the 90% confidence range. More 
inclined NMAGIC models (i — 45°) within the same confi- 
dence range are more radially anisotropic, i.e. /3 ~ 0.6. The 
degree of orbita l an isotropy is consisten t with the previous 
analysis of |N09| and IDeason et all l|2012h . 



6.3 The wider context: intermediate-luminosity 
ellipticals 

NGC 4494 belongs to the group of intermediate-luminosity 
elliptical galaxies with rapidly falling velocity dispersion pro- 
file (NGC 8 21, NGC 3379, NGC 4494 NGC 4697), dubbed 
"naked" by IRomanowskv et al.l (|2003l ) because of the un- 
usually low dark matter content reve aled by their favoured 
dynamical models. The analysis of ICoccato et al.l (2009) 
identified a larger group of early-type galaxies (NGC 821, 
NGC 3377, NGC 3379, NGC 449 4, NGC 4564, NGC 4697) 
with strongly decreasing v ms = \/a 2 ± v 2 profiles. 

Detailed dynamical models have now been ob- 
tained for all of t h e gala xies in the original sample of 
IRomanowskv et al.l |2003l ). including observational con- 
straints from the highe r ord er mo ments of the LOSVD 
(see the summary in | N09l and IWeiimans et al . 2009; 
iForestell fc GebhardtfcOld . for NGC 821). On the whole, the 
results suggest that these intermediate-luminosity galaxies 
are inconsistent with the previous claims of little to no dark 
matter halo. 

We now compare the results that we derived for 
NGC 449 4 in th is w ork, w ith the findings previously ob- 
tained by IDL08I and |DL09| for NGC 4697 and NGC 3379 
using the same modelling technique. Indeed, the photome- 
try and kinematics of these three galaxies are not very dif- 
ferent. They all have intermediate values of luminosity and 
stellar mass, similarly low values of the central velocity dis- 
persion ~ 150 — 210k ms~ 1 , and similarly falling v Iuls profiles 
l|Coccato et al.ll2009h . 

Fig. [23] shows the comparison of the NMAGIC models 
obtained for NGC 4697, NGC 3379, and NGC 4494. For 
NGC 4494, we plot the range of best-fitting edge-on mod- 
els as determined in Section [S] For the other intermediate- 
luminosity ellipticals, no rigorous A\ 2 analysis was per- 
formed, but a range of valid models was determined based 
on the likelihood of PNe. 

As can be seen in the bottom panel of Fig. 1231 the 
orbital structure of the three intermediate luminosity el- 
lipticals is radially biased, in agreement with previous 



(e.g. Merritt & Ohl 1997 


; iRix et al. 19971; Gerhard et all 


19981: van der Marel et all 


1998FlMatthias & Gerhard 1999; 


Gebhardt et al. 2000; Saelia et al.l l2000h. with predictions 


of the monolithic collapse 


scenario (van Albada 1982). and 



with simulations o f both binary- mergers (e.g. Gerhard Il98ll ; 
iDekel et al.l 120051; iThomas et al.l l2009a[) and cosm ological 
mergers (e.g. lAbadi et al.l2006l ; IOnorbe et al.l2007l ). though 
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Figure 23. From top to bottom: as a function of radius, circular 
velocity normalized by its value at IR C , dark matter fraction, and 
anisotropy parameter of the range of valid NMAGIC models ob- 
tained fitting the data of NGC 44 94 (thi s work, violet), NGC 4697 
(IDL08I . orange), and NGC 3379 <DL09t green). 



the remnants of bi nary mergers can exhibit a variety of or- 
bital distributions IjNaab et al.ll2006h . 

NGC 4494 is characterized by a milder level of radial 
anisotropy than the other two galaxies, consistent with the 
fact that i ts projected v e locity dispersion is higher (see 
Fig. 15 of ICoccato et al.l l|2009l )) and its surface density 
profile steeper. Recent simulations have shown that dissi- 
pational processes in wet mergers may decrease the level 
of rad ial anisotropy (e.g. iNaab et al.ll2006l ; iThomas et al.l 
l2009bh , and also explain many observed features such 



nents, and extra light at small radii (e.g. 


Mihos & Hcrnquist 


1994 


: ISDringel et al.l 120051; ICox et all 


2006; Jesseit et al. 


2007 


; Hookins et al.ll200ct iHoffman et all 201011. 



All three galaxies have been classified as fast rotators 
|Emsellem et al.ll201ll). For NGC 4494 its c entral kinemat- 
ically decoupled core dKrainovic et aLlboill). and the small 
value of Ar at large radii ( Coccato et alj 20091 ) . suggest a 
complicated rotation structure. Such t ransitions in Ar to- 
wards large radii (f or other cases, see iProctor et all 2009; 



ICoccato et al. 20091) may be a signature of merger events 



(e.g. IHoffman et al.llioTa V The preference for edge-on mod- 



els found in Section 15.31 also suggests that only a small frac- 
tion of the kinetic energy in this galaxy is in rotation. 

The top panel of Fig. [23] shows the circular velocity 
curves normalized at 1R C . At large radii, the behaviour of 
the circular velocity is not so different, partially due to the 
assumption of logarithmic halo profiles. The greatest dis- 
crepancies can be seen in the central regions, and are proba- 
bly due to different imprints left by baryonic processes dur- 
ing galaxy formation. 

This is also supported by the fact that the dark mat- 
ter fractions of these three galaxies are remarkably differ- 
ent (middle panel of Fig. 123 p , while their global circular ve- 
locity curves are more similar. The dark matter fraction of 
NCG 4494 has much higher values, particularly inside 3R C , 
than those of NGC 4697 and NGC 3379. Such a high dark 
matter fraction could b e a consequence of a merge r event 
(like those advocated bv lProctor et al.l (|2009h or lFllI ). which 



might have scattered dark matter into the in ner regions of 
the galaxy (|Oser et al.ll2012l ; iHilz et al.ll2012l ). 



7 CONCLUSIONS 

We presented dynamical models for the intermediate lumi- 
nosity elliptical galaxy NGC 4494, fitting photometric and 
kinematic observational data to investigate its mass distri- 
bution and orbital structure. Our extended kinematic data 
include d the recently ava ilable integrated light spectra in 
slitlets (|Foster et alj 201 ll ) and hun dreds of planetary nebu - 
lae velocities reaching out to ~ 7R e (|Napolitano et al.l2009h . 
We used the x 2 -made-to-measure particle code NMAGIC to 
construct axisymmetric models for various inclinations ex- 
ploring a large sequence of gravitational potentials. 

In parallel, we carried out a parameter estimation study, 
investigating how well the characteristic parameters of dark 
matter halos can be recovered via NMAGIC modelling of the 
available observational data. For this, we used Monte Carlo 
simulations of NGC 4494-like mock galaxies to determine 
the confidence regions around the best-fitting model. These 
confidence bands were then used to discriminate the range 
of valid models for NGC 4494. 

Our main results can be summarized as follows: 



• Given the observational data of NGC 4494 and 
our NMAGIC modelling set-up, Monte Carlo simulations 
showed that the usual la (70%) confidence level corresponds 
to a relative difference in the merit function AG = 26 
(Ax 2 = 59). At 90% confidence level, AG = 32 and 
Ax 2 ~ 70. These differences are much larger than the com- 
monly used A\ 2 values based on varying a small number of 
model parameters. 

• Under the assumption that the dark matter halo has 
a logarithmic parameterization, the best-fitting models for 
NGC 4494 determined within these confidence levels have 
an approximately flat total circular velocity curve outside 
~ 0.5-R e , with v c (3R e ) ~ 220kms~ 1 . The larger variation 
is in the dark matter circular velocity curve, rather than 
in the stellar one. The dark matter fraction of the mod- 
els within 70% confidence level is about 0.6 ±0.1 at 5R e , 
and they are embedded in concentrated dark halos (ro ~ 1- 
2R C ) with circular velocity ~ [160-230]kms _1 . With this 
large dark matter fraction, the stellar mass-to-light ratio is 
consistent with the value predicted by iDeason et all (|2012h 
from stellar population models with a Chabrier IMF. The 
discrepan cy with the more diffuse halo found by the Jeans 
models of iNapolitano et all l|2009T ) is likely due to the addi- 
tional slitlets constraints, which suggest a milder drop of the 
velocity dispersion than shown by the long-slit kinematics. 

• The edge-on models provide the best fits to the avail- 
able observational data, but the inferred dark halo param- 
eters (ro,vo) do not depend sensitively on the assumed in- 
clination. We explored a sequence of flattened dark matter 
halos, which however did not change our main result: edge- 
on models provide better fits to the data than models with 
lower inclinations. 

• The orbital anisotropy of the stars is increasingly ra- 
dial from the center outwards, but the amount of radial 
anisotropy is smaller than what was found in similar previ- 
ous works for the two other intermediate luminosity ellipti- 
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cals with rapidly falling, quasi-Keplerian velocity dispersion 
profiles, NGC 3379 and NGC 4697. 

• Comparing the halos of all three intermediate- 
luminosity, quasi-Keplerian ellipticals modelled with the 
same made-to-measure particle technique (NGC 3379, 
NGC 4697 and NGC 4494), we conclude that they have 
similar global circular velocity curves and outer dark mat- 
ter halos. NGC 4494 shows a particularly high dark matter 
fraction inside ~ 3R e , and a strong central concentration 
of baryons. These differences are probably related to the 
detailed interplay between baryons and dark matter in the 
processes which shaped these galaxies. 
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